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1  Summary 

The  objective  of  the  project  was  to  improve  efficiency  of  accurate  meth¬ 
ods  for  predicting  radiation  characteristics  of  antennas  mounted  on  large 
platforms. 

The  computational  challenge  in  this  type  of  problems  consists  mainly  of 
the  necessity  of  a  detailed  modeling  of  the  antenna  and,  at  the  same  time, 
describing  the  much  larger  platform  and  its  interaction  with  the  antenna 
system.  While  even  a  complex  antenna  subsystem  can  be  solved  without 
much  difficulty  by  the  available  rigorous  MoM  methods,  computation  of  the 
antenna  interaction  with  the  platform  may  be  out  of  reach  of  the  present 
fast  solvers,  and,  at  the  same  time,  asymptotic  high-frequency  solution  tech¬ 
niques  may  not  be  sufficiently  versatile  and  accurate  for  application  to  re¬ 
alistic  problems.  In  addition,  in  problems  of  antenna  radiation  it  would  be 
highly  desirable  to  develop  efficient  direct  solution  methods,  which  would 
allow  predicting  the  radiation  distribution  in  the  entire  angular  range. 

In  this  context,  the  goal  of  the  present  project  was  to  make  progress 
in  developing  methods  alternative  to  the  presently  prevailing  fast  iterative 
solution  techniques.  Some  approaches  along  these  lines  include  Refs.  [1,  2, 
3,  4,  5,  6]. 

Our  work  was  concentrated  on  an  attempt  of  devising  a  viable  computa¬ 
tional  scheme  based  on  utilization  of  numerically  constructed  basis  functions 
defined  on  large  supports,  and  characterized  by  strongly  collimated  radiation 
patterns;  in  the  following  we  refer  to  them  as  “directional  basis  functions” 
(DBFs).  The  concept  of  the  solution  scheme  was  to  model  the  antenna  using 
the  conventional  MoM  techniques,  to  parameterize  the  platform  in  terms  of 
DBFs,  and  to  describe  interactions  partly  iteratively  (for  interactions  corre¬ 
sponding  to  multiple-scattering  “bounces”),  and  partly  by  means  of  direct 
solution  methods  (for  the  remaining  interactions,  including  propagation  of 
creeping  waves). 

In  the  present  project,  Monopole  Research  concentrated  on  development 
of  the  DBF  construction  techniques,  while  the  Yale  University  pursued  the 
area  of  research  related  to  direct  solution  methods. 

The  main  paxt  of  our  work  was  to  implement  the  concept  of  DBFs  as 
a  mathematically  precise  description  of  well  collimated  beams  in  terms  of 
specific  current  distributions  generating  such  radiation  patterns.  This  ap¬ 
proach  may  be  regarded  as  a  rigorous  realization  of  the  concept  of  “rays” 
in  high-frequency  asymptotics  -  a  mechanism  not  taken  into  account  in  the 
conventional  impedance  matrix  compression  methods.  It  is  hoped  that  such 
an  approach  will  result  in  a  highly  sparse  impedance  matrix  (dominated  by 
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couplings  due  to  matching  collimated  beams  radiated  and  received  by  the 
DBFs),  which,  in  addition,  might  be  amenable  to  the  application  of  direct 
solution  methods. 

We  stress  that  it  appears  difficult  to  achieve  the  above  goals  by  means 
of  constructing  analytically  parameterized  DBFs  based  on  asymptotic  high- 
frequency  methods,  for  at  least  two  reasons:  (a)  the  complexity  and  com¬ 
putational  cost  of  the  algorithms  grows  rapidly  with  the  order  of  the  high- 
frequency  scattering  mechanisms,  and  (b)  in  general,  for  a  given  support  size, 
the  resulting  basis  functions  are  not  guaranteed  to  generate  radiated  beams 
of  optimal  angular  concentration;  in  fact,  radiation  patterns  generated  by 
such  basis  functions  tend  to  have  high  “side-lobes”  (unless  the  behavior  of 
the  function  near  the  boundary  is  regularized  by  using  smooth  partition  of 
unity  or  similar,  rather  complex,  techniques). 

In  this  context,  the  principal  advantages  of  the  numerical  approach  to 
the  construction  of  directional  basis  functions  presented  here  is  its  relative 
independence  of  the  scatterer  geometry  and  the  optimal  angular  collimation 
of  the  radiation  patterns  of  the  DBFs. 

The  main  results  we  obtained  in  this  effort  are  as  follows: 

•  We  have  been  able  to  implement  a  numerical  scheme  for  construction 
of  DBFs,  applicable  to  realistic  surface  geometries,  and  providing,  in 
practice,  near-optimal  angular  collimation  of  radiation  patterns.  We 
confirmed  that  the  angular  widths  of  the  pattern  scale  in  the  expected 
way  with  the  support  size.  We  illustrate  the  result  of  the  algorithm 
on  a  number  of  cases  involving  realistic  surface  geometries  (see  Sec¬ 
tion  2.3). 

•  We  found  the  numerical  cost  of  DBF  construction  is  be  relatively  high, 
and,  as  expected,  rapidly  growing  with  the  size  of  the  basis  function 
support.  However,  on  the  basis  of  our  effort  under  a  parallel  contract, 
where  we  worked  on  a  multilevel  scheme  of  DBF  construction,  we 
expect  the  cost  can  be  reduced  to  the  level  making  the  procedure 
practical  in  realistic  problems  (see  Section  2.4  for  a  brief  description). 

•  We  analyzed  a  number  of  scattering  problems  discretized  by  means 
of  the  constructed  DBFs.  We  encountered  here  difficulties  related  to 
ill-conditioning  of  the  impedance  matrix  in  the  DBF  space  resulting, 
essentially,  from  the  insufficient  linear  independence  of  the  DBFs.  The 
conditioning  problems  seem  to  be  of  rather  fundamental  nature,  and 
can  be  ascribed  to  an  inherent  conflict  between  conditioning  and  an¬ 
gular  collimation,  which  can  be  summarized  in  a  statement  that 
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a  complete  set  of  collimated  “band-limited”  beams  (i.e.,  beams 
generated  by  a  spatially  limited  source)  must  involve  angularly 
overlapping  beams,  and  the  overlap  may  lead  to  ill-conditioning. 

One  can  also  view  construction  of  DBFs  as  a  problem  closely  analogous 
to  antenna  (or  antenna  array)  synthesis,  i.e.,  an  inherently  ill-posed 
inverse  scattering  problem. 

In  view  of  these  difficulties,  we  directed  a  large  part  of  the  effort  to 
the  analysis  of  the  conditioning  problem: 

•  We  carried  out  a  detailed  analysis  of  the  origin  of  the  ill-conditioning 
problems,  and  initiated  efforts  at  their  resolution.  These  efforts  are 
discussed  in  Section  2.5: 

—  In  Section  2.5.2  we  describe  an  attempt  of  modifying  the  opera¬ 
tors  whose  eigenfunctions  ( “radiation  modes” )  form  a  basis  from 
which  the  DBFs  are  constructed.  The  modification  is  expected 
to  allow  us  to  limit  the  range  of  eigenvalues  of  the  operators, 
and  thus  to  improve  conditioning  of  the  DBF  system.  An  inter¬ 
esting  feature  of  this  approach  is  that  the  operators  in  question 
describe  spatial  concentration  of  the  “radiation  modes” ,  and  al¬ 
low  construction  of  basis  functions  concentrated  in  both  angular 
(essentially,  Fourier)  and  configuration  space.  However,  this  tech¬ 
nique  requires  that  the  supports  of  DBFs  significantly  overlap  in 
space,  which  may  again  deteriorate  conditioning.  We  have  not 
yet  reached  a  definite  conclusion  on  the  viability  of  the  approach. 

—  In  Section  2.5.3  we  report  on  an  implementation  of  a  simple  itera¬ 
tive  solution  scheme,  in  which  the  DBF-space  impedance  matrix 
is  used  merely  as  a  compressed  matrix  representation.  In  this 
approach  the  “near-field”  part  of  the  matrix  is  not  discretized  in 
terms  of  DBFs,  which  results  in  a  significantly  better  condition¬ 
ing. 

•  In  the  area  of  the  development  of  fast  direct  solvers,  the  Yale  University 
group  designed  and  implemented  a  fast  direct  solver  for  objects  in 
both  two  and  three  dimensions  that  are  long  and  thin.  An  important 
feature  of  the  scheme  is  its  effectiveness  in  both  high  and  low-frequency 
environment. 

One  of  the  principal  tools  in  the  development  of  algorithms  of  this  type 
is  the  concept  of  “skeletonization” .  It  was  shown  that  under  certain 
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conditions  (highly  relevant  to  the  design  of  scattering  algorithms) ,  in¬ 
troducing  a  randomized  element  into  skeletonization  schemes  leads  to 
radically  improved  efficiency  and  reliability. 

In  summary,  we  believe  that  our  STTR  Phase  I  effort  contributed  to  the 
development  of 

•  a  well  defined  mathematical  procedure  for  specification  and  numerical 
construction  of  directional  basis  functions, 

•  new  efficient  algorithms  (with  randomized  elements)  applicable  to  fast 
direct  solution  methods. 

We  anticipate  that,  as  soon  as  the  encountered  difficulties  are  resolved,  the 
above  developments  will  form  important  building  blocks  for  an  efficient  ad¬ 
vanced  antenna  pattern  prediction  software. 

2  Numerical  construction  of  directional  basis 
functions  (DBFs)  (Monopole  Research  report) 

2.1  General  features  of  the  approach 

We  consider  a  problem  of  electromagnetic  radiation  or  scattering  on  a  per¬ 
fectly  conducting  (closed  or  open)  surface  S.  In  this  case  the  source  of 
radiation  are  the  surface  vector  (electric)  currents,  tangential  to  the  surface; 
we  denote  these  sources  s(r),  r  €  <S. 

The  problems  of  acoustics  and  some  problems  of  elastodynamics  can 
be  formulated  by  simply  replacing  the  vector  sources  s(r)  by  scalar  ones, 
s(r).  More  general  problems  of  elastodynamics  will  require  using  vector 
displacement  fields. 

In  order  to  specify  the  problem  of  constructing  DBFs,  we  construct  a 
set  of  bounded  connected  areas  (“patches”)  covering  the  surface  <S  (we 
may  have  to  allow  for  some  overlap  between  the  patches). 

Our  goal  is  to  construct,  for  each  patch  A4j,  a  set  of  basis  functions 
\ka)(r),  supported  on  the  patch,  which  would  generate  a  complete  set  of 
maximally  collimated  radiation  patterns.  The  two  requirements  which  would 
ensure  an  efficient  numerical  solution  scheme  are: 

(i)  a  high  angular  collimation,  and  the  resulting  high  impedance  matrix 
sparsity,  and 
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(ii)  a  high  degree  of  linear  independence  of  basis  functions,  and  thus  a 
reasonable  conditioning  of  the  resulting  basis  and  of  the  impedance 
matrix  in  that  basis. 

As  we  indicated  above,  these  two  requirements  are,  to  some  extent,  in  mutual 
conflict:  a  finite  patch  size  implies  that  the  generated  radiation  patterns 
are  band- limited  in  the  angle  (or  in  Fourier  space);  therefore,  the  angular 
distributions  must  have  some  overlapping  “tails”,  and  cannot  be  strictly 
linearly  independent. 

In  the  context  of  high-frequency  scattering  it  is  useful  to  separate  the 
possible  current  distributions  into  a  set  of  band- limited  functions  (oscillating 
no  more  rapidly  than  allowed  by  the  wave  number  k),  and  the  remaining 
ones,  involving  higher  Fourier  transforms,  limited  only  by  discretization. 
The  class  of  band-limited  functions  is  sufficient  to  correctly  reproduce  far- 
field  interactions,  and  it  is  only  this  set  which  we  attempt  to  represent  in 
terms  of  DBFs. 

In  the  following  we  give  a  brief  overview  of  our  approach  to  construction 
and  utilization  of  DBFs,  present  the  main  results,  and  describe  attempts  at 
resolving  the  encountered  difficulties. 

2.2  Overview  of  the  solution  procedure  with  DBFs 

In  order  to  establish  the  main  concepts,  we  give  here  a  short  description  of 
our  approach  to  numerical  construction  of  the  DBFs  and  their  use  in  solving 
radiation  and  scattering  problems. 

2.2.1  DBF  construction 

We  consider  a  frequency-domain  radiation  or  scattering  problem  in  free 
space,  characterized  by  the  wave  number  k  and  the  wavelength  A  =  2n/k. 
The  asymptotic  (transverse)  field  radiated  in  the  direction  q  by  the 
considered  tangential  vector  sources  s  on  a  given  patch  M.  =  is 

F00(q)=B(q)  J d2re-ifc^'rs(r)  ,  (2.1) 

M 

where 


n(q)  =  /  —  q  q  ■  (2.2) 

projects  the  field  onto  the  plane  orthogonal  to  the  radiation  direction.  (In 
the  electromagnetic  scattering  problem  with  a  perfectly  conducting  surface 
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the  source  s  is  the  surface  electric  current,  and  is  the  electric  field.)  The 
total  radiated  power  is  then 

P  =  |d2g|F0O(q)|2  (2.3) 

S2 

where  S 2  denotes  the  unit  two-sphere  of  directions  q,  and  the  angular  inte¬ 
gral  is  normalized  to  unity,  f  d2q  =  1. 

Similarly,  we  define  the  power  radiated  into  a  certain  solid  angle1  fl  C  S2 
as 

m  =  /  d^lF^q)!2.  (2.4) 

n 

In  terms  of  the  sources  s  the  power  (2.4)  can  be  expressed  as 

P{Q)  =  j d2rdVs(r)  •  fin(r  -  r')s(r')  =  (s,  HM{9)  s)  ,  (2.5) 

M 

where  the  r.h.s.  defines  the  integral  operator  HM(Cl)  and  the  inner  product 
(• ,  •),  and  where  the  dyadic  integral  kernel  hn  is  given  by 

Mr)  =  J  d2<7  (/  —  qq  ■)  elfcq  r  .  (2.6) 

n 

In  the  case  of  the  total  power  (0  =  S 2)  we  have 

P  =  J  d2rd2r's(r)  • /i(r  -  r')  s(r')  =  (s,  HMs)  ,  (2.7) 

M 

with 

h(r)  =  J  d2q(I -qq-)eikq  r  =  ~[l  - ^{kr)]  +  j2(kr)rr- ,  (2.8) 

S2 

where  the  last  expression  (with  the  spherical  Bessel  functions  j„)  is  related  to 
the  imaginary  part  of  the  Green  function  of  the  usual  electric-field  equation. 
In  the  following  we  refer  to  the  operators  and  HM  as  the  angular 

1An  important  theoretical  and  practical  question  is  how  to  specify  a  reasonable  set  of 
such  angles.  We  discuss  this  problem  in  Section  2.3,3. 
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and  total  power  operators  for  the  given  surface  patch  M.  We  assume  in  the 
following  that  the  operators  act  in  the  space  L2(M)  of  square  integrable2 
tangential-vector  functions  s(r)  on  M. 

We  now  formulate  the  radiation  collimation  problem  in  terms  of  the 
Rayleigh  quotient 


rn  (s.ffjutn)  s) 


(2.9) 


representing  the  fraction  of  the  power  radiated  into  the  solid  angle  Q  relative 
to  the  total  power.  By  construction, 


0  <  fo[s]  <  1  .  (2.10) 

The  source  distribution  s  generating  a  radiation  optimally  collimated 
within  the  given  angle  0  is  the  function  s  maximizing  the  ratio  £n[s].  We 
construct  a  set  of  DBFs  associated  with  the  given  surface  patch  M  simply 
as  a  set  of  source  (current)  distributions  sa,  a  =  1, 2, . . .  ,  of  sufficiently  high 
angular  concentration, 


1-T<fn[s«]<l,  (2.H) 

where  0  <  r  <£.  1  is  an  appropriately  defined  “power  leak  tolerance”  (we 
refer  to  the  (small)  quantity  1  —  £n[sa]  as  the  power  leak  -  the  fraction  of 
the  power  radiated  outside  the  angle  O). 

Clearly,  the  number  of  solutions  corresponding  to  angularly  concentrated 
radiation  patterns  will  be  limited  by  the  size  of  the  patch,  and  will  also 
depend  on  its  shape.  For  sufficiently  regular  patches  it  will  be  of  order  of 
the  Shannon  number  of  the  patch  (proportional  to  its  area  in  the  units  of 
the  wavelength  squared),  i.e.,  to  the  number  of  linearly  independent  band- 
limited  functions  supported  on  the  patch.  It  follows  that  the  set  of  DBFs 
(understood  as  generating  collimated  radiation  distributions)  is  not  sufficient 
to  represent  arbitrary  source  distributions  on  the  patch  (limited  only  by 
discretization).  The  complete  set  of  basis  functions  must  also  include  current 
distributions  weakly  radiating  and  involving  high  Fourier  transforms.  For 
simplicity  of  notation  we  will  denote  all  these  functions  by  'FQ,  understanding 
that  only  a  part  of  them  is  associated  with  collimated  radiation  patterns. 

2The  square  integrability  assumption  may  be  too  weak  in  the  case  of  electromagnetics, 
as  is  it  knowm  that  requirements  of  locally  finite  energy  of  the  fields  usually  restrict  currents 
to  smaller  spaces,  typically  Sobolev  spaces.  This  question  requires  further  analysis  in  the 
present  context. 
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Maximization  of  the  angular  concentration  factor  (2.9)  leads  (by  means 
of  functional  differentiation  with  respect  to  s)  to  the  stationary  point  con¬ 
dition,  which  has  the  form  of  the  generalized  eigenequation 

J  d V  ha( r  -  r')  sa(r/)  =  L  j  dV  h{ r  -  r')  sQ(r')  ,  (2.12a) 

M  M 

or 

HM(n)sa  =  ZaHMsa,  (2.12b) 

where  a  labels  the  eigenvalues  and  eigenfunctions.  In  order  to  find  the 
source  distribution  of  the  best  collimated  radiation  pattern,  we  have  to  find 
the  eigensolution  to  the  highest  eigenvalue  £. 

The  angular  collimation  problem  just  defined  belongs  to  the  category  of 
“concentration  problems” ,  which  have  been  investigated  by  Slepian,  Landau, 
and  Poliak  [7,  8,  9,  10]  for  band-limited  functions  on  the  real  line,  and  later 
generalized  to  functions  of  several  variables  [11,  12,  13,  14,  15,  16,  17]. 

As  follows  from  the  definitions  of  the  kernels  (2.6)  and  (2.8),  the  operator 
HM(Cl)  (and  thus  HM)  is 

1.  self-adjoint, 


HM(Q)  =  ,  (2.13) 

since 

Mr~r')  =  ftn(r'-r)  ,  (2.14) 

2.  positive  semi-definite  (from  Eqs.  (2.4)  and  (2.5)),  and 

3.  Hilbert-Schmidt  (hence  compact ),  since 

Tr  {HXt(Sl)HM(Q))  =  j  d2r  dV  |ftn(r  -  r')|2  <  oo  .  (2.15) 

M 

The  above  properties  of  the  operators  ensure  that  their  spectra  are  discrete, 
non-negative,  of  finite  degeneracy  for  positive  eigenvalues,  and  with  the  only 
possible  condensation  point  at  zero. 

A  conventional  approach  to  solving  the  generalized  eigenproblem  (2.12) 
is  to  first  solve  the  eigenproblem  defined  by  the  r.h.s.  operator  HM,  i.e., 

^ M  ^  u  • 


(2.16) 
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The  spectrum  is  non-negative  and  discrete,  and  the  eigenfunctions  can  be 
chosen  orthonormal, 

(f^  .  U  =  J dVf^r)  ■  f„(r)  =  <^„  .  (2.17) 

M 

Physically,  the  eigenfunctions  are  (square-integrable)  current  distribu¬ 
tions  on  the  patch  M.,  to  which  we  refer  as  “radiation  modes”.  The  corre¬ 
sponding  eigenvalues  are  the  powers  radiated  by  these  currents, 

=  (2.18) 

(cf.  Eq.(2.7));  we  assume  in  the  following  that  these  eigenvalues  are  indexed 
in  the  descending  order. 

By  expanding  the  solutions  sa  of  Eq.(2.12)  in  terms  of  the  eigenfunctions 

Mr)  =  Evyr)’  (2-19) 

we  obtain  the  discrete  (matrix)  eigenequation 


^2  ’ 
V 

^v)  Visa  —  ^  v  /xi/  Vva  ^aV^y^ia  * 

V 

(2.20) 

By  defining 

y»a  —  \J y iia, 

(2.21a) 

and 

—  \/y~  y/v~  ’ 

(2.21b) 

Eq.(2.20)  can  be  brought  to  the  conventional  matrix  eigenequation  form 


^2  y^a  ~  ZaViM  (2.22) 

V 

with  a  positive  semi-definite  Hermitian  matrix  [/i(f2)^„],  and  with  the  re¬ 
sulting  orthonormal  eigenvectors, 

Va  Vb  —  ^  /  V /ia  y»b  ~~  ^ ab 


(2.23) 
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(the  superscript  H  denotes  Hermitian  conjugation) .  In  terms  of  these  eigen¬ 
vectors,  the  eigensolutions  of  the  generalized  eigenequation  (2.12),  i.e.,  the 
DBFs,  are  then,  finally, 


’M1’)  =  sa(r)  =  ^2  Vf/x(r)  =  2 

/i  n 


^af/i(r)  ■ 


(2.24) 


The  appearance  of  the  factors  1  /y/vji  in  this  expression  is  of  crucial 
importance,  as  it  may  easily  cause  ill-conditioning  of  the  system  of  eigen¬ 
functions  sa.  In  particular,  since 

=  (2.25) 

u  ^ 


the  DBFs  vI'a  would  have  been  orthonormal ,  were  it  not  for  the  factors 
1/VV  In  their  presence,  and  if  the  eigenvalues  span  a  large  range,  the 
conditioning  of  the  set  of  eigenfunctions  {’F0}  is  expected  to  have  a  larger 
condition  number. 

In  fact,  since  the  operator  HM  is  infinitely  dimensional  and  compact,  its 
spectrum  must  have  a  concentration  point  at  rj  =  0.  Therefore,  unless  we  can 
cut  off  the  sum  over  the  eigenvalues  at  a  reasonable  threshold,  7/a  >  7?mjn,  the 
system  of  DBFs  may  be  seriously  ill-conditioned.  We  expect  similar  features 
to  persist  also  for  the  discretized,  finite-dimensional  operators.  Therefore, 
we  devoted  a  large  part  of  work  in  this  project  to 


•  analyzing  in  more  detail  the  origin  of  ill-conditioning  of  the  DBF  sys¬ 
tem,  and 

•  investigating  possible  alternative  formulations  of  the  problem  (in  par¬ 
ticular,  modifications  of  the  relevant  operators)  with  the  aim  of  im¬ 
proving  the  conditioning. 

These  questions  are  discussed  in  the  following  subsections,  after  we  an¬ 
alyze  the  role  of  conditioning  of  the  DBF  system  in  solving  scattering  prob¬ 
lems. 


2.2.2  Solution  of  radiation  and  scattering  problems  in  terms  of 
DBFs 

For  definiteness  and  simplicity  we  consider  electromagnetic  scattering  prob¬ 
lem  described  by  the  electric-field  integral  equation  (EFIE)3  on  the  consid- 

3However,  the  general  solution  scheme  discussed  here  does  not  depend  on  the  specific 
properties  of  the  integral  operator  and  the  equation. 
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ered  perfectly  conducting  surface  S. 

j  dV  G( r  -  r')  J(r')  =  -Ein(r)  ,  for  r  G  «S  ,  (2.26) 

s 


where  Em(r)  is  the  incident  field,  and  G  is  the  Green  function 


G(r)=s(r)-pvff(r)V- 

with  the  Helmholtz  equation  Green  function 

eik  |r| 


0(r)  = 


4  7r  |r | 


(2.27) 


(2.28) 


In  order  to  simplify  the  presentation,  we  do  not  consider  the  conventional 
MoM  discretization  of  the  above  equations,  but  rather  directly  compare 
discretizations  in  terms  of  the  sets  of  the  orthonormal  modes  and  in 
terms  of  DBFs  \E,a. 

In  the  first  case  we  represent  the  discretized  equations  in  the  matrix  form 


Ax  =  b 


(2.29) 


where  b  represents  the  incident  field,  x  the  solution,  and  A  the  impedance 
matrix,  all  having  block  structure  associated  with  the  patches  of  the  surface. 
If  Galerkin  discretization  is  used,  the  current  on  the  patch  Mi  is  expanded 
as 


J(0(r)  =  £  if  fW(r)  for  r  €  Mt  ,  (2.30) 

elements  of  the  vector  b  are  projections  of  the  incident  field  on  the  basis 
functions  on  the  patch, 

bf  =  -<f«  ,  Ein)  =  -J  d2r  fG)(r)  •  Ein(r)  ,  (2.31) 

and  the  elements  of  the  matrix  A  blocks  are 

a<il>  =  (tf\GfP)=  j  d2r  j  dV  tf  (r)  •  G(r  -  r')  f <#)  (r')  .  (2.32) 

Mi  Mj 
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For  reasonably  regular  patch  geometries  and  typical  MoM  discretizations 
with  spatial  resolution  of  order  10  points  per  wavelength,  condition  numbers 
of  the  matrix  blocks  are  moderate,  say,  of  order  103  for  patch  size  of  about 
4  A. 

The  discretized  equations  in  the  DBF  space  are 

Ax  =  b  (2.33) 

with  the  impedance  matrix  elements 

=  <*«  ,  GV^)  =  J  d2r  J  dV^fr)  •  G{ r  -  r')  ^(r')  ,  (2.34) 

M{  Mj 

and  expressions  for  the  current  and  the  r.h.s.  analogous  to  Eqs.  (2.30)  and 
(2.31). 

It  follows  from  Eq.(2.24)  that  the  DBF  matrix  blocks  are  transforms  of 
the  matrix  blocks  in  the  modes  of  the  power  operator, 

®aP  =  53  yiP  >  (2.35a) 

flU 

or 

A(ij)  =  Y®  H  AM  Y ^  ,  (2.35b) 

where  Y ®  and  Y^  are  transformation  matrices  consisting  of  the  elements 
y W  and  for  patches  M;  and  Mj.  In  terms  of  the  transformation  matrix 
Y 

Y  =  diag[y^  •  •  •  yW]  (2.36) 

(where  p  is  the  number  of  patches)  we  have  then 

A  =  YhAY,  (2.37a) 

b  =  Ynb ,  (2.37b) 

x  =  Yx,  (2.37c) 


i.e.,  we  can,  formally,  transform  the  original  equation  (2.29)  into  Eq.(2.33), 
solve  Eq.(2.33)  (with  a  sparse  matrix  A),  and  transform  the  solution  x  to  x. 
However,  since  the  transformation  matrix  blocks  may  have  large  condition 
numbers  (due  to  factors  1 /  y/vii  in  Eq.(2.24)),  the  transformed  matrix  A 
blocks  (in  particular,  the  diagonal  ones)  may  also  be  ill-conditioned,  and 
solution  of  the  system  (2.33)  may  be  practically  impossible. 

In  Section  2.5  we  address  the  conditioning  problem  in  mode  detail  and 
describe  our  attempts  at  its  resolution. 
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2.3  Construction  of  DBFs  -  numerical  examples 

2.3.1  Construction  of  eigenfunctions  of  the  power  operator 

To  illustrate  construction  of  eigenfunctions  of  the  power  operator  for 
patches,  we  consider  a  curved  surface  patch  of  sizes  3  A  x  3  A  (Fig.  1)  located 
approximately  in  the  (x,  y)  plane,  and  discretized  with  n  =  1384  MoM 
unknowns  (associated  with  the  edges).  The  singular  values  a ^  (i.e.,  square- 
roots  of  the  eigenvalues  rj)  of  the  power  operator  are  shown  in  Fig.  2  (in 
this  particular  example  we  used  the  singular-value  decomposition  algorithm 
instead  of  the  equivalent  eigenvalue  analysis) . 

10° 

10"’ 

10"* 

10° 

S  10”4 
to-4 
to-4 
10-T 
1(T* 

0  200  400  600  800  1000  1200  1400 

M 

Figure  2:  Distribution  of  the  sin¬ 
gular  values  of  the  power  operator 
Hm  of  the  patch  of  Fig.  1. 

Some  of  the  “radiating  modes”  (eigenfunctions  of  H j^)  are  shown  in 
Figs.  3-6.  Figs.  3  and  4  show  the  current  distribution  and  the  radiation 
pattern  for  the  eigenfunction  associated  with  the  highest  eigenvalues  r]1 ,  and 
Figs.  5  and  6  the  same  quantities  for  a  very  weakly  radiating  mode  with  the 
eigenvalue  7j400  ~  10-10  r)l.  The  radiation  patterns  are  computed  at  Gauss- 
Legendre  quadrature  points  on  the  unit  sphere  S2  (for  the  quadrature  order 
L  =  27),  and  plotted  in  an  approximate  Mercator  projection:  the  indices 
0  <  <  2L  label  the  quadrature  points  in  the  range  0°  <  <p  <  360°,  and  the 

indices  0  <  n6  <  L  label  the  quadrature  points  in  the  range  180°  >  9  >  0° 
(tiq  =  0  corresponds  to  the  vicinity  of  the  south  pole,  9  =  180°). 

Generally,  the  current  distributions  and  radiation  patterns  of  the 
strongly  radiating  modes  are  more  regular  than  for  the  weakly  radiating 
ones.  A  more  careful  analysis  shows,  actually,  that  the  strongly  radiating 
modes  are  of  the  type  of  standing  waves.  However,  their  radiation  patterns 


Figure  1:  The  surface  patch  M 
(n  =  1384  MoM  unknowns)  used 
in  the  analysis. 
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are  practically  never  collimated  to  any  appreciable  degree.  For  example, 
although  the  radiation  pattern  of  Fig.  4  is  peaked  in  two  approximately  op¬ 
posite  directions  (6  ~  90°,  4>  ~  90°,  270°),  the  level  of  radiation  outside  the 
peaks  is  not  much  reduced. 

A  characteristic  feature  of  the  weakly  radiating  modes  is  a  large  content 
of  high  Fourier  transforms ,  typically  higher  than  the  incident  wavelength, 
and,  near  the  end  of  the  spectrum,  approaching  the  highest  oscillation  rate 
allowed  by  the  spatial  sampling  (discretization  of  the  surface) .  This  property 
is  consistent  with  the  fact  the  eigenfunctions  associated  with  small  eigenval¬ 
ues  are  not  band-limited. 


Figure  3:  The  real  part  of  the  Figure  4:  The  radiation  pattern  fj 

y-component  of  the  eigenfunction  generated  by  the  eigenfunction  of 

fi,  corresponding  to  the  highest  Fig.  3. 
eigenvalue. 
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#q:  1457/1458 


Figure  5:  The  real  part  of  the  Figure  6:  The  radiation  pattern 

y-component  of  the  eigenfunction  f400  generated  by  the  eigenfunc- 

f4oo,  corresponding  to  a  small  tion  of  Fig.  5. 
eigenvalue. 

2.3.2  Construction  of  DBFs 

As  we  described  in  Section  2.2.1,  construction  of  DBFs  amounts  to  finding 
angularly  well  concentrated  solutions  of  the  generalized  eigenequation  (2.12), 

i.e.,  solutions  corresponding  to  eigenvalues  close  to  unity. 

For  typical  geometries,  the  spectrum  of  the  eigenvalues  £a  consists,  qual¬ 
itatively,  of  three  subsets: 

1.  eigenvalues  close  to  1,  i.e.,  solutions  concentrated  on  f l, 

2.  eigenvalues  close  to  0,  i.e.,  solutions  excluded  from  Cl,  and 

3.  the  remaining  “transition  eigenvalues”,  associated  with  solutions  nei¬ 
ther  concentrated  on,  nor  excluded  from  0. 

The  number  of  solutions  concentrated  on  a  given  region  Cl  grows  with  the 
size  of  the  region.  This  behavior  is  similar  to  that  observed  in  the  case  of 
the  functions  concentrated  on  spherical  cups  [13,  16,  17];  however,  in  our 
problem  it  depends  on  the  subsystem  geometry  and  on  the  shape  of  the 
desired  angular  region. 
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The  minimum  angular  half-width  of  the  region  (in  the  case  the  DBF  is 
concentrated  in  one  region  only)  is  about  20N,  where 

*  “  (2-38) 

is  the  “Nyquist  angle”  for  a  radiating  system  contained  in  a  box  of  side  D  (in 
our  example  0N  ~  11°).  Again,  however,  such  a  statement  is  approximately 
valid  only  for  a  geometry  filling  more-or-less  uniformly  the  entire  box. 

An  important  feature  of  the  geometry-dependent  radiation  patterns  is 
that  in  some  cases,  if  a  high  concentration  is  to  be  achieved,  the  region 
Q  must  consists  of  several  disjoint  areas.  For  example,  currents  on  a  flat 
surface  radiate  symmetrically  with  respect  to  the  surface  plane,  and  thus 
the  region  Q  must  also  be  symmetric.  A  similar  case  is  a  thin,  wire-like, 
object,  for  which  the  admissible  radiation  patterns  must  be  symmetric  with 
respect  to  the  object  axis. 

We  describe  more  details  of  our  solution  procedure  in  the  following  Sec¬ 
tion.  Here  we  only  mention  that  we  start  with  a  set  of  tentatively  defined 
angular  regions  0  and  gradually  expand  each  region  until  the  at  least  one 
eigen-solution  is  sufficiently  well  concentrated,  i.e., 

!-*«<£,  (2-39) 

where  e  is  the  power  leak  tolerance,  related  to  the  criteria  of  accuracy  for 
matrix  elements  and  for  sparsity  of  the  DBF-space  impedance  matrix;  we 
typically  take  e  in  the  range  10-5  to  10-4. 

We  show  below  some  examples  of  DBFs  constructed  for  the  patch  M  of 
Fig.  1,  following  the  procedure  outlined  in  Section  2.2.1. 

In  Fig.  7  we  plot  distributions  of  the  eigenvalues  £a  for  a  set  of  six  regions 
each  consisting  of  two  areas,  concentrated  near  the  poles  ( 0  ~  0°  and 
0  ~  180°).  These  regions  are  created  during  the  DBF  construction  procedure 
in  which  the  initial  region  is  being  gradually  expanded  until  the  condition 
(2.39)  is  met  for  at  least  some  eigen-solutions. 
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Figure  7:  Distribution  of  the  con¬ 
centration  factors  £a  for  a  set  of 
six  angular  regions  Q  of  increas¬ 
ing  sizes,  with  the  indicated  num¬ 
bers  of  directions  within  the  re¬ 
gion  (the  total  number  of  direc¬ 
tions  is  2  L2  =  1458). 


Figure  8:  Distribution  of  the 
power  leaks  1  —  corresponding 
to  the  same  set  of  angular  regions 
as  in  Fig.  7. 


In  the  considered  DBF  construction  our  code  has  stopped  expanding  the 
area  0  after  six  steps,  and  selected  14  DBFs  with  power  leaks  below  the 
assumed  tolerance  e  =  10-5.  The  radiation  patterns  of  the  first  and  last  of 
these  DBFs  are  shown  in  Figs.  9  and  10. 
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#q:  704/ 1458  #q:  755/ 1458 


Figure  9:  The  radiation  pattern  Figure  10:  The  radiation  pattern 
|'I'i|  of  the  first  selected  DBF  for  | ] 4 1  of  the  last  selected  DBF  for 
the  largest  angular  region  Cl  of  the  largest  angular  region  Cl  of 
Fig.  7.  Fig.  7. 

Figs.  7  and  8  show  how  the  number  of  well  concentrated  solutions  in¬ 
creases  with  the  growing  size  of  the  region  Cl.  Typically,  the  minimum 
angular  half-width  of  the  region  is  about  2  0N,  where  0N  =  A/(\/3  D)  is  the 
“Nyquist  angle”  for  a  radiating  system  contained  in  a  box  of  side  D.  In  our 
case  0N  ~  11°. 

However,  in  our  example  we  selected  the  tentative  angular  region  in 
the  direction  of  the  positive  2-axis  (the  north  pole) .  Because  of  the  relative 
flatness  of  the  considered  patch  (Fig.  1)  it  appears  impossible  to  collimate  the 
radiation  in  one  direction  only  -  here  in  the  north-pole  direction  -  without 
radiating  fields  in  the  symmetric  direction.  This  effect  is  seen  in  Fig.  9  and 
especially  in  Fig.  10,  where  there  is  a  significant  radiation  emitted  in  the 
direction  of  the  south  pole  ( n6  ~  0).  Consequently,  the  selected  angular 
region  Cl  is  rather  wide,  and,  finally,  contains  about  one- half  of  all  directions 
(this  is  partly  due  to  the  fact  that  the  density  of  the  quadrature  points  near 
the  poles  is  about  twice  as  high  as  near  the  equator). 

Figs.  11  and  12  show  similar  results  for  a  more  favorable  case  of  the 
tentative  radiation  direction  along  the  rc-axis  ( 9  —  90°,  <p  —  0°).  In  this 
case  it  was  possible  to  obtain  a  much  narrower  radiation  pattern  of  the 
constructed  DBFs.  Our  algorithm  expanded  the  initial  angular  region  in 
four  steps,  and  selected,  eventually,  four  solutions  satisfying  the  condition 
(2.39). 
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Figure  11:  Distribution  of  the  con¬ 
centration  factors  (a  for  a  set  of 
four  angular  regions  Cl  of  increas¬ 
ing  sizes,  with  the  indicated  num¬ 
bers  of  directions  within  the  re¬ 
gion. 

ftq:  126/ 1458 

B0001_00056_00001  - 

0.0001  - 

0.001 


Figure  12:  Distribution  of  the 
power  leaks  1  —  £Q  corresponding 
to  the  same  set  of  angular  regions 
as  in  Fig.  11. 
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Figure  13:  The  radiation  pat¬ 
tern  |^i  |  of  the  first  selected  DBF 
for  the  last  angular  region  Cl  of 
Fig.  11. 


Figure  14:  The  radiation  pat¬ 
tern  |  ^4 1  of  the  last  selected  DBF 
for  the  last  angular  region  Cl  of 
Fig.  11. 


2.3.3  Determination  of  angular  concentration  regions  for  DBFs 

In  the  previous  Section  we  showed  results  of  constructing  DBFs  with  radia¬ 
tion  patterns  concentrated  in  prescribed  angular  regions  Cl.  We  now  briefly 
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describe  the  procedure  we  use  to  find  such  regions  in  a  way  consistent  with 
the  geometry  of  the  patch,  and  with  the  resulting  physical  restrictions  on 
the  admissible  radiation  patterns. 

Our  goal  is  to  determine  a  given  number  AT  of  DBFs  \Ha  concentrated 
with  the  tolerance  e  (Eq.(2.39))  on  as  small  as  possible  angular  regions  fiQ 
(since  the  size  of  the  regions  0  controls  the  sparsity  of  the  transformed 
impedance  matrix  blocks) .  The  number  A f  of  DBFs  must  be  at  least  equal 
to  the  Shannon  number  of  the  patch,  i.e.,  the  approximate  number  of  linearly 
independent  band-limited  functions  (current  distributions)  supported  on  the 
patch;  this  number  is 


Af(M)  =  2n 


IM 

A2  ’ 


(2.40) 


where  \M\  is  the  area  of  the  patch,  and  the  additional  factor  2  comes  from 
the  two  polarizations  of  the  electromagnetic  field. 

Carrying  out  the  task  of  minimizing  the  angular  regions  in  a  rigorous  way 
would  lead  to  an  extremely  complex  combinatorial  problem,  which  would 
be,  in  practice,  impossible  to  solve.  We  resort,  therefore,  to  approximate, 
heuristically  motivated  (and  based  on  numerical  experience)  methods. 

We  start  with  specifying  a  set  of  approximately  uniformly  distributed 
directions  and  associated  regions  Cl  in  which  we  would  like  the  DBFs  to 
radiate.  We  take  the  number  of  regions  as 


N0  = 


pAf 


n0 


(2.41) 


where  p  1  is  a  “redundancy  factor”  and  ne  is  the  desired  average  number  of 
well-concentrated  eigen-solutions  per  region.  We  adjust  this  number  based 
on  the  considerations  of  numerical  efficiency,  since  the  cost  of  constructing 
more  eigensolutions  for  a  single  region  is  lower  than  constructing  them  for 
separate  regions;  on  the  other  hand,  specifying  too  large  regions  may  deteri¬ 
orate  angular  concentration.  We  found  that  a  reasonable  number  is  of  order 
10. 

The  redundancy  factor  accounts  for  the  likely  possibility  that  (e.g.,  due 
to  symmetries  of  the  geometry)  not  all  created  DBFs  will  be  sufficiently 
linearly  independent;  it  also  compensates  for  the  unlikely  cases  where  the 
DBF  construction  fails  (i.e.,  no  concentrated  solutions  are  obtained  for  any 
reasonably  small  region  0  -  see  below).  The  number  of  solutions  ne  depends 
on  the  size  of  the  region,  which  implies  that  we  have  to  choose  “tentative” 
region  sizes  accordingly.  An  approximate  value  of  the  opening  half-angle  9n 
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of  the  region  fl  is 

0n~2v/^*N  =  ^£^.  (2.42) 

Having  determined  the  set  of  “tentative”  angular  regions  Cl  of  angular 
widths  0^,  we  construct  DBFs  for  each  region  Cl  independently.  As  the  first 
approximation,  we  evaluate  the  spherical-cap  concentration  function  [13], 
say  x(q)i  of  the  assumed  angular  width  concentrated  on  Cl.  We  then 
project  that  angular  distribution  on  the  space  spanned  by  radiating  modes 
u^.  As  a  result,  we  obtain  a  modified  angular  distribution  x(q)  whose  shape 
exhibits  the  features  of  the  geometry;  e.g.,  for  a  geometry  with  a  symmetry 
plane,  the  modified  distribution  will  have  two  peaks,  symmetric  with  respect 
to  the  plane.  We  use  the  projected  angular  distribution  to  determine  the 
updated  region  Cl  as  the  set  of  directions  in  which  the  values  of  x  exceed 
a  given  threshold.  We  then  solve  the  eigen-problem  for  the  region  Cl ,  and, 
if  the  required  concentration  conditions  are  not  met,  we  keep  modifying 
(usually  expanding)  the  region  Cl  until  we  obtain  at  least  ne  acceptable 
solutions.  The  region  Cl  is  being  updated  based  on  the  behavior  of  the  best 
concentrated  solutions:  we  include  in  Cl  directions  in  which  the  obtained 
radiation  patterns  exceed  a  threshold  value. 

The  above  procedure  of  finding  the  regions  Cl  usually  terminates  after 
few  (typically  about  5)  steps,  depending  on  the  selected  direction,  initial  size 
of  the  region,  thresholds,  etc.  In  rare  cases  the  region  Cl  expands  until  it 
covers  the  full  solid  angle  (or  most  of  it) .  These  “failed”  solutions  are  likely 
to  be  eliminated  in  the  stage  of  selecting  the  optimal  subset  of  linearly 
independent  DBFs  (as  described  below). 

In  our  example  we  constructed  Na  =  134  initial  angular  regions,  and  we 
obtained,  on  the  average,  5  acceptable  solutions  per  region,  i.e.,  the  total  of 
710  DBFs  -  about  1.3  times  more  than  the  required  number.  There  were  no 
“failed”  construction  cases,  and  the  power  leaks  for  the  DBFs  ranged  from 
1  — £  =  3.54e— 8  to  1  —  £  =  9.99e-6,  with  the  average  1— £  =  3.79e— 6.  The 
angular  regions  covered  from  about  8  %  to  53  %  of  the  quadrature  points, 
with  the  average  of  18  %. 

2.4  A  multilevel  scheme  of  DBF  construction 

One  of  the  practical  difficulties  in  numerical  DBF  construction  is  its  compu¬ 
tational  cost,  rapidly  increasing  with  the  patch  size  and  the  number  of  MoM 
unknowns  n  on  the  patch.  The  main  contribution  to  the  cost,  growing  as 
n3,  is  due  to  the  eigenvalue  (or  singular-value)  analysis.  At  the  same  time, 
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sufficiently  large  patch  sizes  are  necessary  in  order  to  achieve  an  impedance 
compression  which  would  would  be  competitive  with  those  of  more  conven¬ 
tional  matrix  compression  methods. 

In  order  to  alleviate  this  difficulty,  we  have  investigated  (within  a  dif¬ 
ferent  contract)  a  possibility  of  a  hierarchical,  multi-level  scheme  of  DBF 
construction,  in  which  DBFs  on  large  patches  are  expressed  as  superposi¬ 
tions  of  DBFs  on  smaller  sub-patches.  Here  we  briefly  summarize  the  main 
results,  in  order  to  indicate  that  construction  of  DBFs  on  large  patches  (of 
sizes  of  tens  of  wavelengths)  may  be  feasible. 

Our  approach  to  the  multi-level  DBF  construction  is  analogous  to 
antenna-array  synthesis,  in  which  the  array  radiation  pattern  is  a  superpo¬ 
sition  of  radiation  patterns  of  its  elements  (DBFs  of  the  sub-patches)  with 
coefficients  ensuring  an  optimal  angular  collimation  of  the  array  pattern. 
Such  an  optimization  problem  results,  again,  in  a  generalized  eigenequa- 
tion  for  the  coefficients  multiplying  the  sub-patch  DBFs.  In  a  more  detailed 
analysis  of  the  problem  we  found  that  a  nearly  optimally  collimated  DBF  ra¬ 
diation  into  a  certain  solid  angle  Q  can  be  constructed  by  taking  into  account 
only  those  DBFs  for  sub- patches  which  radiate  in  directions  overlapping  the 
angular  region  0;  this  fact  significantly  lowers  the  computational  cost  of  the 
construction. 

As  a  numerical  example,  we  present  the  DBF  computation  for  a  set  of 
four  adjacent  level- 1  patches  of  sizes  about  3 A  x  3 A  of  a  curved  surface, 
forming  a  single  level-2  patch  of  size  about  6A  x  6A. 

Typical  results  are  given  in  the  Table  I  below.  The  third  column  gives  the 
percentage  of  the  solid  angle  covered,  on  the  average,  by  the  DBF  radiation 
pattern,  defined  by  requiring  that  the  fraction  of  power  radiated  outside  the 
DBF  angular  region  is  less  than  10  .  The  total  matrix  compression  is  the 

product  of  the  factor  in  the  fourth  column  and  an  additional  compression 
(ranging  here  from  0.40  to  0.10)  due  to  elimination  of  a  part  of  weakly 
radiating  modes. 
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Table  I 

DBF  collimation  and  construction  costs 


patch  size 

unknowns 

|n|/4ir 

compression 

time/unknown  time/unknown 

single  level 

multi-level 

1.5  A 

330 

45% 

0.80 

0.2  s 

- 

3.0  A 

1300 

18% 

0.14 

2.0  s 

2.0  s 

6.0  A 

5200 

6% 

0.01 

16.0  s 

2.2  s 

The  Table  shows  scaling  of  the  width  of  the  radiation  pattern  of  the  DBF 
with  its  support  size,  the  resulting  sparsity  of  the  far-field  impedance  matrix 
blocks,  and  is  consistent  with  the  computational  cost  of  DBF  construction 
(per  unknown)  approaches  a  constant,  as  expected  on  theoretical  grounds 
(the  total  cost  of  constructing  all  DBFs  should  scale  with  the  number  of 
unknowns  N  as  0{N  log  N)).  The  times  given  in  the  Table  were  obtained 
in  a  computation  on  a  single  AMD  Athlon  processor. 

Figs.  15  -  18  below  visualize  the  behavior  of  a  level-2  DBF  constructed 
using  all  level- 1  DBFs,  and  only  the  subset  of  the  angularly  overlapping 
DBFs. 


#q:  129/5832 


#q:  129/5832 


h_0001_MQ  -  h  0001_MQ  - 

0.0001  -  0.0001  - 

0.001  -  0.001  - 


Figure  15:  The  radiation  pattern 
of  a  level-2  DBF  constructed  using 
all  level-1  DBFs. 


Figure  16:  The  radiation  pattern 
of  Fig.  15  plotted  in  logarithmic 
scale. 
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#q:  134/5832  #q:  134/5832 

h_0001_MQ  - —  h_0001_MQ 

0.0001  - -  0.0001 


Figure  17:  The  radiation  pattern  Figure  18:  The  radiation  pattern 
of  a  level-2  DBF  constructed  using  of  Fig.  17  plotted  in  logarithmic 
only  a  subset  of  overlapping  level-  scale. 

1  DBFs. 

2.5  Analysis  of  the  conditioning  problem 

We  give  here  an  account  of  our  analysis  of  the  ill-conditioning  difficulties 
associated  with  the  DBFs,  and  of  our  attempts  at  their  resolution.  We  start 
with  indicating  a  relation  between  band-limited  source  distributions  and 
far-field  couplings;  we  then  proceed  to  discussing  modifications  of  the  power 
operators  which  lead  to  improved  spectra  of  band-limited  eigenfunctions, 
and  to  describing  an  alternative  iterative  solution  scheme  which,  to  large 
extent,  circumvents  the  conditioning  problem. 

2.5.1  Computation  of  “far-field”  matrix  blocks  in  the  DBF  space 

In  Section  2.2.2  we  expressed  the  matrix  elements  of  the  impedance  ma¬ 
trix  in  the  DBF  space  A  in  terms  of  the  matrix  elements  in  the  basis  of 
the  “radiation  modes”.  Similarly,  the  matrix  A  can  be  evaluated  by  first 
computing  matrix  blocks  in  the  MoM  basis,  and  then  transforming  them  to 
the  DBF  space.  However,  such  a  procedure  is  computationally  expensive 
(the  cost  per  transformed  blocks  is  0(n3),  where  n  is  the  average  number 
of  unknowns  associated  with  the  patch) . 

For  this  reason,  in  the  implementation  of  the  algorithm  we  used  an  alter¬ 
native  expression  for  the  DBF  matrix  elements,  based  on  the  Fast  Multipole 
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Method  (FMM),  which  allows  us  to  compute  the  blocks  of  the  matrix  A 
for  sufficiently  distant  patches  at  a  much  lower  cost,  proportional  to  n2  for 
patches  of  sizes  larger  than  few  wavelengths.  We  mention  this  algorithm  at 
this  point,  since  it  also  allows  us  to  find  upper  bounds  on  the  error  in  the  far- 
field  matrix  elements  due  to  possible  approximations  done  in  construction 
of  DBFs. 

We  consider  DBFs  4/^  and  4^  associated  with  two  patches,  Mi  and 
Mj,  such  that  the  distance  between  the  centers  c2  and  c-  of  the  smallest 
spheres  enclosing  the  patches  is  larger  than  the  sum  of  the  radii  R(  and  Rj 
of  the  spheres,  i.e., 


cij  =  \cij\  =  lci  -  cjl  >  Ri  +  Rj  ■  (2-43) 

In  this  case  we  say  that  the  patches  are,  mutually,  in  their  “far-field”  range. 
We  can  then  utilize  the  well-known  Fast  Multipole  Method  (FMM)  expres¬ 
sion  for  the  matrix  element  (2.34), 

=  J  d2Q^a\k^)c.  •T(cii,q)n(q)$Jj)(A;q)c.  ,  (2.44) 

52 

where  II(q)  is  given  by  Eq.(2.1), 

¥p{kq)Ci  =  J  d2re-ik^r-cJyM{r)  (2.45) 

is  the  Fourier  transform  of  the  DBF  4'^  relative  to  the  center  Cj  of  the 
patch  (and  similarly  for  the  other  DBF),  and  T  is  the  FMM  diagonal 
form, 


• »  ^max 

T(c ij,  q)  =  ^  ^  (2i  +  1)  ie  h{1}  (k  Cij)  P<(cy  •  q)  .  (2.46) 

e=o 

Here  h^  is  the  spherical  Hankel  function  of  the  first  kind,  is  the  Legendre 
polynomial,  and  the  truncation  in  the  angular-momentum  sum  is  approxi¬ 
mately  Lmax  >k(Ri  +  Rj). 

The  essential  feature  of  the  representation  (2.44)  is  that  the  matrix  ele¬ 
ment  is  expressed  in  terms  of  the  Fourier  transforms  of  the  basis  functions 
with  the  arguments  not  larger  than  k.  This  means  that  the  far-field  ma¬ 
trix  elements  are  sensitive  only  to  the  band-limited  components  of  the  basis 
functions. 
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The  above  property  of  the  representation  (2.44)  allows  us  to  obtain  an 
upper  bound  on  the  error  in  the  DBF  matrix  elements  due  to  an  approxima¬ 
tion  in  the  DBF:  if  the  original  basis  function  4/^  is  modified  to  +  A'l'a  \ 

then,  to  the  first  order  in  A’Fa  ,  the  variation  in  the  matrix  element  is 
bounded  by  the  expression 

lAaifl2  <  M(|Cii|)P[A4/«]  P[^}  ,  (2.47) 

where  M(  |cf  •  |)  is  a  function  depending  only  on  the  distance  between  the 
patches,  and  the  factors  P[---]  are  total  powers  radiated  by  the  current 
distributions  A 'J'q ^  and  defined  as  in  Eqs.  (2.1),  (2.2),  and  (2.3). 
Eq.(2.47)  allows  us  to  relate  the  error  in  the  matrix  elements  to  the  eigen¬ 
values  r/;(  of  the  radiation  modes  included  and  neglected  in  the  computation 
of  DBFs. 

2.5.2  Alternative  power  operators  and  their  spectral  properties 

As  we  discussed  in  Sections  2.2.1  and  2.2.2,  the  spectra  and  the  eigenfunc¬ 
tions  of  the  power  operator  HM  are  the  decisive  factors  controlling  the 
conditioning  of  the  constructed  set  of  DBFs. 

We  have  investigated  possibilities  of  improving  the  conditioning  by  tak¬ 
ing  advantage  of  some  arbitrariness  in  the  definition  of  the  angular  concen¬ 
tration  of  the  radiation  patterns.  This  arbitrariness  is  due  to  the  fact  that, 
if  a  radiation  pattern  is  weighted  by  a  smooth  function,  it  remains  concen¬ 
trated  independently  of  the  shape  of  that  function  (provided  its  variation  is 
small  over  the  concentration  region). 

We  found  that,  by  choosing  appropriate  weighting  functions  it  is  possible 
to  construct  such  modified  angular  power  distribution  that  the  correspond¬ 
ing  power  operators  (defined  as  in  Section  2.2.1)  have  properties  similar  to 
those  of  the  “concentration  problem”  operators,  such  as  giving  rise  to  prolate 
spheroidal  wave  functions  and  their  generalizations,  and  which  we  encoun¬ 
tered  in  the  angular  concentration  problem.  More  specifically,  the  modified 
power  operators  define  concentration  of  the  radiation  modes  in  space  (on 
the  patch  surface). 

The  difference  in  the  spectra  of  the  original  and  the  modified  power 
operators  is  shown  in  Fig.  19  for  a  curved  patch  M  of  size  3A  x  3A.  We  show 
the  spectrum  of  the  original  operator  HM  (normalized  to  =  1),  and  two 
modified  operators  denoted  HMl  and  HM2. 
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Figure  19:  Spectra  of  the  original  power  operator  Hm  and  modified  opera¬ 
tors  Hmi  and  Hm2  for  a  patch  of  size  3A  x  3A. 

The  spectra  of  the  modified  operators  are  similar  to  those  of  the  angular 
concentration  problem  (Figs.  7  and  11):  the  highest  eigenvalues  are  close  to  1 
(corresponding  to  current  distributions  well  concentrated  on  the  patch),  and, 
after  a  transition  region,  the  eigenvalues  rapidly  (exponentially)  approach 
0,  corresponding  to  spatially  “deconcentrated”  distributions. 

By  using  the  modified  power  operators  instead  of  the  original  ones,  we 
are  able  to  limit  the  range  of  the  eigenvalues  (since  there  are  now  more 
eigenvalues  close  to  1  in  the  relevant  part  of  the  spectrum);  this  fact  improves 
conditioning  of  the  DBF  set.  On  the  other  hand,  the  spatially  concentrated 
eigenfunctions,  being  band-limited,  have  to  smoothly  vanish  at  the  patch 
boundaries,  which  requires  utilizing  overlapping  patches,  and  causes  deteri¬ 
oration  of  conditioning.  The  eventual  result  depends  on  the  choice  of  the 
amount  of  overlap,  and  on  the  procedure  of  eliminating  a  part  of  the  lin¬ 
early  dependent  DBFs;  we  have  not  yet  been  able  to  find  appropriate  criteria 
ensuring  reasonable  DBF  conditioning. 

2.5.3  An  alternative  scattering  problem  solution  scheme 

In  Section  2.2.2  we  formulated  the  scattering  problem  by  discretizing  the 
integral-equation,  alternatively,  in  terms  of  the  “radiation  modes”  and  in 
terms  of  DBFs.  In  close  analogy,  the  conventional  MoM  matrix  equation 


Ax  =  b 


(2.48) 
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can  be  transformed  to  the  DBF  form 

Ax  =  b  (2.49) 

(Eq.(2.33))  with 

A  =  YhAY,  (2.50a) 

b  =  Ye  b  ,  (2.50b) 

x  =  Yx  ,  (2.50c) 

where  Y  is  the  transformation  matrix  expressing  DBFs  as  linear  combina¬ 
tions  of  the  MoM  basis  functions.1 2 * 4 

However,  ill-conditioning  of  the  system  of  DBFs  (and  thus  the  trans¬ 
formation  matrix  F),  causes  ill-conditioning  of  the  impedance  matrix  A, 
which  may,  in  practice,  preclude  solving  the  problem  directly  (by  matrix 
inversion).  In  fact,  in  the  problems  we  were  attempting  to  solve,  we  were 
not  able  to  obtain  stable  solutions  by  means  of  the  LU  decomposition  of  the 
matrix,  and  even  the  much  more  expensive  Singular  Value  Decomposition 
(SVD)  algorithm  was  yielding  ambiguous  results. 

In  view  of  these  difficulties,  we  implemented  an  alternative,  iterative  so¬ 
lution  scheme,  in  which  the  DBFs  are  utilized  only  to  construct  a  compressed 
representation  of  the  impedance  matrix.  The  procedure  is  as  follows: 

1.  Separate  the  MoM  impedance  matrix  into  near-field  and  far-held 
blocks  (associated  with  pairs  of  patches),  defined  according  to  the 
FMM  criteria  (Section  2.5.1), 

A  =  An  +  Af  .  (2.51) 

In  our  implementation  we  define  patches  by  simply  partitioning  the 
geometry  into  cubic  boxes.  The  near-held  matrix  block  are  then  those 
coupling  nearest-neighbor  boxes,  and  the 

2.  Transform  the  far-held  part  of  the  matrix  to  the  DBF  space, 

Af  =  Y*A{Y  (2.52) 

(in  the  practical  implementation  we  compute  the  blocks  of  A  directly 
rather  than  by  transformation,  as  mentioned  in  the  last  footnote). 

4In  a  practical  implementation  direct  computation  of  the  matrix  A  and  its  transforma¬ 
tion  would  be  too  costly,  and  the  transformed  equation  has  to  be  constructed  in  a  more 
economical  way.  In  our  implementation  we  compute  MoM  matrix  blocks  and  transform 
them  only  in  the  near-field  range,  and  evaluate  the  far-field  matrix  blocks  by  means  of  the 
FMM  algorithm,  as  described  in  Section  2.5.1,  Eq.(2.45). 


3.  Compute  the  inverse  transformation  matrix 


=  Y'1 


(2.53) 


(this  matrix  is  also  block-diagonal). 

3.  Solve  iteratively  the  original  equations  (2.48)  with  the  matrix  A  rep¬ 
resented  as 


A  =  An  +  Rh  A{  R  (2.54) 

and  with  the  matrix- vector  multiplication  y  =  Ax  implemented  as  the 
sequence  of  operations 


y  :=  Rx  , 

(2.55a) 

V'=A tV  » 

(2.55b) 

y~RHy, 

(2.55c) 

y:=Anx  +  y. 

(2.55d) 

Thus,  the  above  scheme  uses  the  DBF  representation  to  form  a  compressed 
representation  (2.54)  of  the  original  matrix  A:  in  Eq.(2.54)  An  is  block- 
sparse,  R  is  block-diagonal,  and  blocks  of  A{  are  sparse,  provided  the  colli- 
mation  of  radiation  patterns  of  the  DBFs  is  sufficient. 

While  the  procedure  avoids  transformation  of  near-field  matrix  blocks, 
there  remains  the  problem  of  inverting  the  poorly  conditioned  transforma¬ 
tion  blocks  Y.  In  practice,  we  found  that  task  manageable,  and  were  able  to 
obtain  reliable  iterative  solutions  consistent  with  the  original  MoM  solutions. 

3  Developments  in  the  formulation  of  fast  direct 
solver  (Yale  University  report) 


Our  work  under  this  contract  can  be  roughly  subdivided  into  two  parts: 
development  of  tools  specific  to  the  problems  to  be  addressed,  and  the  work 
involving  collateral  issues.  As  often  happens,  in  a  number  of  cases  collateral 
issues  had  to  be  addressed  in  order  to  deal  effectively  with  the  principal  sub¬ 
ject  of  research.  Below  is  a  brief  discussion  of  both  types  of  issues  addressed 
by  the  project. 
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In  addition  to  these  investigations,  we  have  been  reporting  results  ob¬ 
tained  recently  (but  prior  to  the  start  of  the  current  contract).  These  can 
be  found  in  [21],  [22],  [23], 

1.  We  designed  and  implemented  a  fast  direct  solver  for  objects  in  both  two 
and  three  dimensions  that  are  long  and  thin.  This  work  involved  the  “fast” 
element  of  the  algorithm,  the  rapidly  convergent  discretizations  in  various 
situations,  and  the  required  asymptotic  analysis.  An  important  feature  of 
the  scheme  is  its  effectiveness  in  both  high  and  low-frequency  environment. 
This  work  has  been  reported  in  [18];  the  scheme  has  been  generalized  to  two 
dimensions,  and  the  algorithm  is  being  implemented. 

2.  One  of  spin-offs  of  this  activity  is  an  observation  that  given  a  set  of  n  very 

general  bounded  functions  , tpn  (with  a  finite  n),  there  exist  sta¬ 

ble  n-point  interpolation  and  quadrature  formulae  exact  on  the  linear  space 
spanned  by  the  functions  ,  y>2 ,  •  •  •  ,  <pn-  This  is  a  somewhat  unexpected 
observation,  and  its  proof  depends  entirely  on  linear  algebraic  arguments; 
we  have  also  implemented  numerical  algorithms  for  the  construction  of  such 
interpolation  and  quadrature  formulae  applicable  in  many  situations  of  prac¬ 
tical  interest  (in  addition  to  the  uses  of  such  schemes  in  numerical  scattering 
theory).  These  results  are  reported  in  [19]. 

3.  One  of  the  principal  tools  in  the  development  of  algorithms  of  this  type  is 
the  concept  of  “skeletonization”  (see,  for  example,  [22]).  Recently,  it  turned 
out  that  under  certain  conditions  (highly  relevant  to  the  design  of  scat¬ 
tering  algorithms),  introducing  a  randomized  element  into  skeletonization 
schemes  leads  to  radically  improved  efficiency  and  reliability.  These  results 
are  reported  in  [20]). 


*  4- 
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